function result = CS1_0_extraction_2D_Phi(data_id , extra)
if nargin <2
    extra = 1;
end
pathname = 'raw data/';
[raw_data, ~] = load_3ds(['Grid Spectroscopy' num2str(data_id,'%.3d') '.3ds'], pathname,nan);
flag=0;
for i=1:length(raw_data)
    if strcmp(raw_data(i).channel,'LockinXchannel')
        dIdV1 = raw_data(i).map;
        bias = raw_data(i).e;
        rx = raw_data(i).rx.*1e9*0.9;
        ry = raw_data(i).ry.*1e9*0.9;
    end
    if strcmp(raw_data(i).channel,'Lockin2Xchannel') || strcmp(raw_data(i).channel,'Input6')
        dIdV2 = raw_data(i).map;
        flag=1;
    end
end
if flag==1
    dIdV = dIdV2;
    d2Id2V = dIdV1;
else
    dIdV = dIdV1;
    d2Id2V = dIdV1.*0;
end

[nr,nc,ne] = size(dIdV);

for ic = 1:nc
    for ir = 1:nr
        %%%% gaussfilt then pick local maximum %%%% 
        [~ , tmp] = max( imboxfilt(dIdV(ir,ic,:),3) );
        phi_1(ir,ic) = bias(tmp);
        n_fit = 2;
        tmp_fit(:) = d2Id2V(ir,ic,:);
        tmp_y = tmp_fit( max(1,tmp-n_fit) : min(tmp+n_fit,ne) );
        tmp_x = bias( max(1,tmp-n_fit) : min(tmp+n_fit,length(tmp_fit)) );
        reg = [ones(length(tmp_x),1) tmp_x'] \ tmp_y';

        phi(ir,ic) = -1 * reg(1)/reg(2);
        if extra == 0 && (tmp < 2 || tmp > ne-1)
            phi(ir,ic) = bias(tmp);
        end
    end
    
end

%find the center of the defect
phi3 = smoothdata( smoothdata(phi_1','gaussian',nr/2)' ,'gaussian',nc/2);
if data_id==100
    phi3 = smoothdata( smoothdata(phi_1','gaussian',20)' ,'gaussian',20);
end

%find the center of the defect: then, mask the center area
mask = zeros(nr,nc);
mask( round(nr*2/5) : round(nr*3/5) , round(nc*2/5) : round(nc*3/5) ) = 1;

[tmp_phi ,indy] = min(phi3.*mask);
[~ , indx] = min(tmp_phi);
indy=indy(indx);

result.phi = phi;                       %in eV
result.rx = rx - rx(1,indx);            %in nm
result.ry = ry - ry(indy,1);
result.lx = max(max(rx))-min(min(rx));  %in nm
result.ly = max(max(ry))-min(min(ry));
result.nx = size(rx,2);
result.ny = size(rx,1);
result.nc = nc;
result.nr = nr;
result.indy = indy;
result.indx= indx;
result.lb = 25.66 / 6^0.5;
result.EC = 56.1 * 6^0.5;

%build distance map
tmpy = meshgrid(1:nr)'-indy;
tmpx = meshgrid(1:nc)-indx;
[~,polXY] = cart2pol(tmpy,tmpx);
%flattern the matrix
polXY2= reshape(polXY,[],1);
phi22 = reshape(phi,[],1);
%sort
[polXY2_sorted,I] = sort(polXY2);
phi22_sorted = phi22(I);

%group data by pixel
id_p = 1;
pixel = 0:1:nr/2*1.1;
n_grouped=length(pixel);

phi_grouped_count = zeros(n_grouped,1);
phi_grouped=cell(1,n_grouped);
for i=1:nr*nc
    if polXY2_sorted(i)>pixel(id_p)
        id_p = id_p+1;
    end
    if id_p>length(pixel)
        break;
    end
    
    phi_grouped_count(id_p) = phi_grouped_count(id_p)+1;
    phi_grouped{id_p} = [phi_grouped{id_p} , phi22_sorted(i)];
end

result.az_r = pixel * (result.lx/(nc-1));


%%%% remove outlier %%%%
%median is more stable against outliers
for i=1:n_grouped
   for j=phi_grouped_count(i):-1:1
       if abs(phi_grouped{i}(j) - median(phi_grouped{i})) >1e-3 
           phi_grouped_count(i)= phi_grouped_count(i)- 1;
           phi_grouped{i}(j) = [];
       end
   end
end

for i=1:n_grouped
    result.az_phi(i) = mean(phi_grouped{i});
    result.az_phi_err(i) = std(phi_grouped{i});
end
result.phi_grouped_count = phi_grouped_count;
result.avg = mean(result.az_phi(end-10:end));

end